Stochastic thermodynamics of chemical reaction networks 
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For chemical reaction networks in a dilute solution described by a master equation, we define en- 
ergy and entropy on a stochastic trajectory and develop a consistent nonequilibrium thermodynamic 
description along a single stochastic trajectory of reaction events. A first-law like energy balance re- 
lates internal energy, applied (chemical) work and dissipated heat for every single reaction. Entropy 
production along a single trajectory involves a sum over changes in the entropy of the network itself 
and the entropy of the medium. The latter is given by the exchanged heat identified through the 
first law. Total entropy production is constrained by an integral fluctuation theorem for networks 
arbitrarily driven by time-dependent rates and a detailed fluctuation theorem for networks in the 
steady state. Further exact relations like a generalized Jarzynski relation and a generalized Clausius 
inequality are discussed. We illustrate these results for a three-species cyclic reaction network which 
exhibits nonequilibrium steady states as well as transitions between different steady states. 

PACS numbers: 05.40.-a, 05.70-a, 82.39.-k, 82.60.-s 



I. INTRODUCTION 

Networks of coupled chemical reactions in a dilute so- 
lution should be described by a chemical master equation 
whenever fluctuations are relevant due to small numbers 
of at least one of the involved speciesi^^^. As essen- 
tial parameters this equation contains the rate constants 
of all possible reactions. The solution of the chemical 
master equation gives the dynamics of the probability to 
find a certain number of molecules of each species at a 
given time for a given initial condition. The stochastic 
character of such a description implies that for each re- 
alization we can talk about the stochastic trajectory of 
the network by recording, in principle, the time at which 
each particular reaction took place with its concomitant 
change of the number of molecules. Even though such a 
description can be applied to networks in thermodynamic 
equilibrium, it is particularly relevant for driven networks 
in open systems where externally imposed boundary con- 
ditions generate fluxes even at a constant (mean) con- 
centration of species. More generally, such a descrip- 
tion holds true even if the rate "constants" become time- 
dependent, e.g., if the concentration of a chemical species 
can be controlled externally in a time-dependent fashion. 

Quite often, such networks arc effectively at constant 
temperature and pressure, i.e., at thermodynamically 
well-defined conditions. This holds true in particular 
for the small networks discussed in cell biology where 
the few relevant molecules are embedded in the aque- 
ous intracellular solution. Such networks describing, e.g., 
gene regulation 8 -!*^, signal transduction (see, e.g^i) or 
molecular motor a 12 i 13 ' 14 i 15 ' 16 , operate typically in non- 
equilibrium generated by (time-dependent) external me- 
chanical or chemical stimuli. Feedback loops can gener- 
ate cyles as a characteristics of non-equilibrium steady 
states. Taking both the fluctuation aspects and the well- 
defined thermodynamic conditions seriously, the question 
arises whether or not it is possible to give a thermody- 
namically consistent description of such networks on the 
level of a single (stochastic) trajectory. 



Thermodynamics on the level of a fluctuating trajec- 
tory in a driven system may look far fetched at first sight. 
Such an approach, however, has been extremly fruitful 
for somewhat different but conceptually related systems 
consisting of a few degrees of freedom, like colloidal par- 
ticles or (bio)polymersi 7 -. The dynamics of these sys- 
tems in aqueous solution under the action of external 
forces generated mechanically by the tip of an atomic 
force microscope or optically by laser tweezers can be 
typically described by a Langevin equation. For such 
a driven dynamics, a thermodynamically consistent de- 
scription has recently been worked out based on two 
essential ingredients which can also be seen as the es- 
sential assumptions of such an approach, (i) A first- 
law like balance relates external work, internal energy 
and exchanged heat dynamically along each fluctuating 
trajector y 18 ! 19 , (ii) The heat dissipated into the envi- 
ronment leads to a well-defined entropy increase of the 
surrounding medium which is assumed to keep locally 
its constant temperature despite the fact that the few 
degrees of freedom of the system (colloidal particle or 
polymer) embedded into this bath are driven to non- 
equilibrium. For such systems, various exact relations 
have been derived theoretically. The most prominent 
examples are the Jarzynski relatio n 20 i 21 ' 22 i 23 which al- 
lows to extract free energy differences from experiment or 
simulations performed under non-equilibrium conditions 
and the fluctuation theore m 24 ' 25 i 26 ' 27 i 28 ' 29 which quanti- 
fies the probability to observe entropy annihilating tra- 
jectories in the steady state. These relations and various 
ramification a 23 i 30 i 31 have been tested in real and com- 
puter experiments both using mechanical stretching of 
biomolecule o 32 i 33 ' 34 i 35 ' 36 i 37 as well as colloidal particles 
driven by laser trap o 38 i 39 ' 40 i 41 . Since the fundamental 
constituents (a stochastic trajectory, a surrounding heat 
bath and some source of non-equilibrium) of these me- 
chanical systems are also present in chemical networks 
as described above, a similar approach should be feasible 
for them as well. The purpose of this paper is to develop 
this correspondence and explore its consequences in some 
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generality formally. 

Fluctuation theorems have indeed been discussed pre- 
viously both for the steady state of non-equilibrium 
networks by Gaspard and coworker a 6 ' 7 ! 42 and for time- 
dependently driven chemical reactions in Ref4^ with- 
out making the connection to a first-law energy balance. 
Qian and co-workers have developed such an energetic 
perspective mainly on the ensemble leve l 44 i 45 ' 46 i 47 ex- 
tending Hill's classical work 4 ^. On the trajectory level, 
unpublished work by Shibata 4 ^ and our previous studies 
on simple models for enzyms^S, point in this direction. 

The present systematic approach towards a stochas- 
tic thermodynamics^ for chemical networks based on a 
master equation rests on two main ingredients. First, we 
apply the first law to a single trajectory, i.e. we analyse 
a single reaction event in terms of chemical work ap- 
plied by the chemiostats, change in internal energy and 
amount of dissipated heat. This allows us to express the 
exchanged heat in each reaction step in terms of the rate 
constants of the master equation. Second, we define en- 
tropy production also along a single trajectory as a sum 
of two contributions. The heat dissipated in each step 
increases the entropy of the surrounding medium. More- 
over, there is an entropy of the network itself. The sum 
of both changes obeys various exact relations from which, 
e.g., the second law follows trivially for the mean entropy 
production. 

The paper is organized as follows. In Section [Til we 
recall the stochastic modelling of chemical reaction net- 
works by a master equation approach. In Section IIII1 
we define energy, work and dissipated heat on a stochas- 
tic trajectory and develop a consistent noncquilibrium 
first law along a single stochastic trajectory of reaction 
events. In Section ITVl we define stochastic entropy along 
a trajectory of the system and total entropy production 
of system and heat bath. In Section [V] we derive a fluc- 
tuation theorem for the total entropy production from 
which a second-law like statement follows directly. We 
then relate this fluctuation theorem to the Jarzynski re- 
lation. The total entropy production can be split up into 
a housekeeping heat, which accounts for the irreversibil- 
ity due to noncquilibrium steady states, and a quantity, 
which accounts for the irreversibility due to transitions 
between steady state a 30 ' 51 . Both quantities obey a gen- 
eralized fluctuation theorem. In SectionfVlJ we apply our 
results to a paradigmatic reaction network. 
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FIG. 1: Coupling of the system with species Xj, j — 
(l,...,Nj) to the N a particle reservoirs for species A& at 
chemical potential fj, a and to a heat bath at constant temper- 
ature T. 



with 1 < p < N p labeling the single (reversible) reactions. 
We distinguish two types of reacting species. The Xj 
molecules ( j = 1 , . . . , Nj ) are those species whose num- 
bers n = (ni, . . . , Hjv,) can, in principle, be measured 
along a chain of reaction events. In practice, these num- 
bers should be small. The A a molecules (a = 1, . . . , N a ) 
correspond to those species whose overall concentrations 
c a are controlled externally by a chemiostat due to a 
(generally) time dependent protocol c a (r). In principle, 
this implies that after a reaction event has taken place, 
the used A a arc "refilled" and the produced ones are 
"extracted" . These chemiostats have chemical potential 



Tin (c a /c ) 



(2) 



where pP a is the chemical potential for standard condi- 
tions at concentration Co and T is the temperature of the 
heat bath to which both type of particles are coupled, 
see Fig. [T] Throughout the paper we set Boltzmann's 
constant to unity which implies that entropy will be di- 
mcnsionlcss. We assume that the reacting species have 
no internal degrees of freedom. However, internal degrees 
of freedom could easily be treated within our approach 
by labeling different internal states as different species 
and defining "reactions" (transitions) between them. 

The stochiometric coefficients r a , pj, s a and qj enter 
the stochiometric matrix V with entries 



II. MODEL 
A. Chemical master equation 

We consider the general reaction network^ 



p p p 

< = €, - P j 



(3) 



and the stochiometric matrix of the reservoir species U 
with entries 



(4) 



N a Nj N a Nj 

£ rg A a + £ tf.Xj ^ Yl < A « + E ^ 

a—l 3 — 1 a—l j — 1 



For the externally controlled concentrations c Q of A a , we 
(1) will use the vector notation c = (ci, . . . , c/v a ) in the fol- 
lowing. We assume a dilute solution of reacting species in 
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a solvent (heat bath) and therefore the transition prob- 
abilities per unit time for the N p reactions ((T|) take the 
text book formic 



tu£(n, c) = Qkl Y[ (c Q w Q ) r ° • J 



w<L (n, c) = flk p _ Y\ (c Q w Q ) s ° • J 



Uj\ 



(5) 



Uj\ 



(n 3 --g?)!.(fi/ Wi )* 



(6) 

where + denotes a forward reaction, — denotes a back- 
ward reaction and f2 is the reaction volume. The bare 
rates fc p _ are the transition probabilities per unit time 
per unit volume per unit concentration (in terms of 1 /w a 
and 1/coj, respectively) of every educt reactant. Note 
that while uj p _(n, c), in principle, can be measured ex- 
perimentally, the bare rates fc p _ depend on the normal- 
izing volumes u a and u)j whose unique definition, as we 
will see in Sect. IIII Al requires a microscopic Hamilto- 
nian. 

The transition probabilities depend only on the current 
state and therefore define a Markov process with a unique 
master equation 

d T p(n,r) = h+(n- v p ,c) -p(n- v p ,t) 
p 

+ w p _ (n + v p , c) • p(n + v p , r)] 



[w+(n,c) -p(n,T) 

p 

+w f L (n,c) -p(n,r)] 



(7) 



governing the time evolution of the probability distribu- 
tion p(n, r) to have rij molecules Xj at time r. Here, we 
have used the vector notation v p = (v p , . . . ,v p N .) for the 
entries of the stochiomctric matrix. 



B. Stationary states 

Solving the time-dependent master equation ([7]) is out 
of scope for most systems. For constant c a , the long-time 
limit of any initial distribution is a stationary state un- 
der quite general conditions^. The master equation for 
the relaxation process (at constant c a ) can be solved an- 
alytically for linear reaction networks^. We first explore 
whether any given network allows for an equilibrium solu- 
tion p cq (n, c) which obeys the detailed balance condition 

p oq (n,c)w p (n,c) =p cq (n + v p ,c)u) p (n + v p ,c). (8) 

Using (0) and ©, we get 



P cq K c) _^Ln ( \< TT 
p° q (n + vP,c) " k p + 11^*^ 11 



(n, + v p )m< 



3 ^jlu/ 



Inserting a multivariate Poissonian 

?> c ) = IInK s ; 



(10) 



with mean into equation §§§ yields the N p equations^ 

k p 



_± 

kl 



n 



(ii) 



These N p equations (jTTJ) arc equivalent to the system of 
linear equations 

fc p 

E^ = m -^+E< m ( c «^) (12) 



with 



€,j = - In (rijUij/ri) 



(13) 



We now have to distinguish three cases concerning the 
solvability of the system (TT^J) which lead to three different 
physical situations. 



Case I : Equations flS\) have a unique solution. 

Generically, this case requires N p — Nj. The unique 
solution of eqs. (fT2|) corresponds to the unique equilib- 
rium distribution, if all states are accessible from every 
initial state through a sequence of reaction events. The 
equilibrium distribution p0[) can then be written as 



p^(n,c)=Afl[—^n/u Jj r 



(14) 



with normalization factor N = Y[j e " J ■ Taken at face 
value, the form Q14p resembles a grandcanonical distribu- 
tion with one particle energies tj and chemical potentials 
fij = 0. As we will show below, consistency requires 
splitting 6j in a one-particle energy Ej and a non-zero 
chemical potential fij. 

However, even if eqs. (|12[) have a unique solution, not 
all states must necessarily be accessible from every initial 
state through a sequence of reactions due to the discrete 
reaction kinetics. The equilibrium distribution then is 
the projection of the Poissonian ([14]) to the subspace of 
accessible states. 



Case II : Equations H12[l have more than one solution. 

The rank of the stochiomctric matrix V then is smaller 
than the number Nj of species Xj. This is usually the 
case if there are less reactions p than species Xj, i.e. 
N p < Nj. Wc then have constraints for the accessible 
(9) state space depending on the initial distribution. The 
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state space {n} is then separated into an infinite number 
of subspaces. If all possible initial states lie in the same 
subspace, the equilibrium distribution is the projection 
of the Poissonian (|14[) to this subspace with the ej given 
by any solution of (fl~2|) . As is shown in App. |Bl it does 
not matter which solution of (fT2|) is chosen. 

Case III : Equations H12V have no solution. 

The chosen c a then create a genuine noncquilibrium 
stationary state (NESS) violating the detailed balance 
condition. This will usually happen, if there are more 
reactions p than species Xj, i.e., N p > Nj. We then have 
essentially different paths leading to the same final state, 
i.e., there are genuine cycles in the network. 

Even if the network has no equilibrium solution for the 
given c Q , there must be certain values c® q for which eqs. 
(fl~2")l have a solution. For we can always assume to isolate 
the system after preparing an initial state with given c a . 
The network will then relax into an equilibrium state 
with concentrations c?, q . 



III. FIRST LAW 

We want to state the first law of thermodynamics along 
a stochastic trajectory n(r) for an arbitrary network ([1]). 
We therefore need the concept of internal energy, chemi- 
cal work and dissipated heat. 

A. Internal energy 

We start with internal energy and consider an equi- 
librium situation with concentrations c° q which are even 
defined for case III, as discussed above. Then, a solu- 
tion of (fT2")) exists for this set of reservoir concentrations 
c° q . We also assume a microscopic Hamiltonian for the 
dilute solution of Xj and A a molecules in the solvent 
from which an (average) energy Ej and E a for the re- 
acting species Xj and A a , respectively, can be derived as 
shown in App. [AJ The chemical potential |[5J) then can 
be written as 

li a = E a +T In (c a u a ) (15) 

where the volumes ui a are chosen such that E a becomes 
the energy of one molecule A a as defined by a micro- 
scopic Hamiltonian, see App. [A] for details. If, for ex- 
ample, the A a could be treated as an ideal gas without 
internal degrees of freedom, this normalization volume 
would become uj a = exp(— 3/2) where A Q is the ther- 
mal de Broglie wavelength of an A a molecule and the 
factor exp(— 3/2) guarantees that the thermal kinetic en- 
ergy is included in the energy E a of an A a molecule^. 

Even though the rewriting from eq. @ to eq. (fT^l) 
looks trivial, it is a crucial step since it provides the rela- 
tion between the energy E a which will enter the first law 



and the chemical potential p a . The very fact that there 
is a concentration c Q = \/uj a where both become numer- 
ically identical does, of course, not imply that these two 
quantities are similar in general. In fact, the volume fac- 
tor uj a as derived in App. IA1 depends on the microscopic 
details. 

The key point for the subsequent analysis is that the 
energies Ej thus identified in an equilibrium situation 
from statistical mechanics are independent from c cq . 
Since in a dilute solution, there is no interaction between 
the reacting species apart from the chemical reactions, 
the energy Ej of one molecule Xj will neither depend on 
the number of molecules of the other species nor on the 
concentrations c a . Hence, the Ej are well defined even 
for networks in a NESS and for networks driven by time- 
dependent c q (t). Likewise, the definition of the chemi- 
cal potential of reservoir species (fTS"|) still holds true in 
such noncquilibrium situations, since the particle reser- 
voirs are assumed to be large. 

A system energy or internal energy 

E(n(r)) = J2E 3 n J (r) (16) 

i 

along a stochastic trajectory can then be defined both 
in equilibrium and under the specified non-equilibrium 
conditions. 



B. Work and heat 

For a formulation of the first law along a stochastic 
trajectory, we next have to calculate the heat flowing 
from the system into the heat bath. Whenever a reaction 
p takes place in forward direction (r = +1) or backward 
direction (r = —1), i.e., whenever n changes, the energy 
of the system changes by 

AE r - p = AE ± ^ P = J2 EjAn r j' p = ±J2 E J v j- ( 17 ) 

3 3 

where AnJ' p = rv^ is the change in the number of Xj 
molecules due to the reaction r, p. The chemical work 
done by the particle reservoirs in such a reaction step is 

KL> = W±^ m = - £ p a An^ = T E ( 18 ) 

a a 

where An^' p = ru p a is the number of A a molecules trans- 
formed in the reaction ±p. Due to energy conservation, 
the heat flowing from the system into the heat bath fol- 
lows as 

Q r - P = K&r* - ^ = -j2 ^a<" - r J2 M- 

a j 

(19) 

Thus, we have identified the first law of thermodynamics 
for a single reaction event as 

K£m - + Q^. (20) 
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For a finite time intervall [0, t] with given concentra- 
tion protocols c q (t), we sum (fTT|) . (TTg]) and ([TO]) over all 
occuring reactions at times r;, where wc will denote 
forward reactions by r; = +1 and backward reactions by 
r\ = —1. We then get the chemical work 

W chem = J2W ri ' pt = -5^^A*a(ca(^))An2-« 



2 a 



the change in internal energy 



(21) 



(22) 



and the dissipated heat 



with 



and 



z a 

/ rfr yZ E i' h 3 +y2^a(c a {T))h a (23) 



2 



S(t - n)riv 



(24) 



(25) 



Since only energies Ej and chemical potentials fi a enter 
the crucial quantities AE, Q and W, this formulation of 
the first law is also valid in nonequilibrium. 

The present identification of work, internal energy and 
heat depends crucially on the fact that we control the 
concentrations c a externally. This means in particu- 
lar that the chemical work has to be spent for "refill- 
ing" the chemiostats after each reaction step. A some- 
what different identification applies if one considers re- 
laxation (without further external interference) from an 
initially prepared nonequilibrium state with concentra- 
tions c a (r = 0) to the corresponding equilibrium with 
concentrations c® q for r — > oo. In this case, the system 
should comprise the Xj and the A a molecules. Then 
we cannot distinguish between the two types of species 
Xj and A a anymore and we have to label all reacting 
species with Xj, (j = 1, . . . , N a + Nj) leaving no chemio- 

stat species A a . The change in the internal energy of the 
system for a single reaction step then is 



AE 



E 1 



i p E n 



(26) 



and since the first law then involves no external chemical 
work anymore, Q = —AE is dissipated as heat. Our 
results (including those of the following sections) hold for 
this special case, too, if they are applied to the system of 
Xj without chemiostat species. 



C. Relation between statistical and kinetic 
approach 

In equilibrium, the chemical potentials \ij must obey 



E<^ 



E 



< [E a + T hi (d?w a )] 
(27) 

since there is no chemical potential difference between 
connected states in equilibrium. We now compare the 
equilibrium distribution (|14[) to the grandcanonical equi- 
librium distribution of a dilute solution, which can always 
be written as (see App. [X]for details) 



P 



J (n,c eq ) =AfY 



("M) T 



-/3(E j n j - f j, j n j ) 



(28) 



with the energy Ej of one Xj molecule and appropri- 
ate normalization volumes ujj. Using the relations for 
the chemical potentials fij ([27)) and the system of lin- 
ear equations for the tj (|12p . we see that thermodynamic 
consistency constrains the ratio of bare rates k p _ / by 



/3j2v?Ej=ln^-pJ2< E " 



(29) 



where (3 = 1/T is the inverse temperature. 

Two remarks are appropriate. First, in general the 
energies Ej cannot be extracted from the kinetics, i.e. 
the chemical master equation, without knowledge of the 
detailed Hamiltonian since the bare rates k p + _ depend 
on the normalizing volumes u> a and ojj which are defined 
only through a microscopic Hamiltonian. Second, if one 
assumes the k p + _ and the u> a , ujj to be given without 
specifying an underlying microscopic Hamiltonian, one 
has to make sure that the energies E a are chosen such 
that eqs. (|29|) have a solution for the Ej. From another 
viewpoint, given energies E a and Ej constrain the possi- 
ble bare rates for a given thcrmodynamically consistent 
reaction network. An illustration of this fact will be given 
in Sect. [VTI below. 



IV. ENTROPY 

We distiguish an entropy change of the system proper 
from a change in the entropy of the medium (heat bath). 
The total entropy production is the sum of both entropy 
changes. 



A. Medium entropy change 

The dissipated heat Q from the first law corresponds 
to a change in entropy of the heat bath 



As r 



T T 
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along a single trajectory. The change in medium entropy 
due to one occuring reaction p in direction r can be ex- 
pressed through rate constants using the definitions of 
energies (f2U)) and chemiostat chemical potentials (| 15[) as 



^Q r > p = -r/3 (v, M + E ^< j 



r In —5- - r u p a In (c Q w Q ) 



In 



<(n,c) ^ (n/wj) rv 'inj\ 



w_ r (n + rvP,c) 



'Ml 



{31) 



In the third line, we have used the transition probabilities 
([5]) and ([6|). Summing (|31| over all occuring reactions 
in direction 77 at times 77 yields the medium entropy 
change 

As m = ^nln^ r -^r^<m(c a a; a ) (32) 
along a single trajectory. 

B. System entropy 

We will now define an entropy of the system. The total 
entropy production of system and heat bath then is the 
relevant quantity for the second law of thermodynamics. 
Again, we will be guided by the equilibrium case first. 
Equilibrium state functions can be obtained from 
using the partition function 



^ = En 



(33) 



and the grandcanonical potential 

J = -TlnZ= -T^ (n/wj-Je-"^-"'). (34) 



Of particular importance is the equilibrium entropy of 
the system 



S°* = -d T J 



= -]>> oq (n) lnp cc <(n)-E ln 
V 3 

This expression differs by a term 



("Mr 



35) 



involving the degeneracy of the state n from the usual 
Shannon entropy 

S Sh = - ^>(n) Inp(n) = (- Inp(n)) (37) 



of a Markov process. Relaxing the equilibrium 
constraint, we define the ensemble entropy in non- 
equilibrium as 

S(t) = -^p(n,r) (lnp(n,r) -hxg n ) 

11 

= (-lnp(n,r)+ln 5n ). (38) 

We will now generalize this ensemble expression to a 
stochastic entropy valid along a single trajectory for both 
equilibrium and noncquilibrium situations following^. 
We define the stochastic system entropy 

S (r) = -lnp(n(T),T)+ S (n(r)), (39) 

where p(n(r),r) is the solution of the master equation 
([7|) taken along the stochastic trajectory n(r) and 



V n i( T ) ! 



(40) 



is the internal entropy of the state n due to the degen- 
eracy. The ensemble entropy (|38[) then is the average of 
the stochastic entropy ([39|) . 



C. Total entropy production 

The total entropy production along a stochastic tra- 
jectory 



As to 



As + As m 
p(n(0),0) 



In 



p(n(t),t) 



E ln 



l-pl 



E r ' ln lir -E ri E< ln ( c « w ") ( 41 ) 



then is the sum of the changes of system entropy ([39"]) and 
medium entropy (|3"2"|) , where As = s(t) — s(0). It will be 
shown in Sect. IV Al below that, in fact, Astot is inde- 
pendent of the normalizing volumes tu a and tUj. Hence, 
this quantity does not require explicitly a microscopic 
Hamiltonian but is rather determined by the full rates 
wi. _ (n, c) entering the master equation. As an aside, we 
note that after averaging the stochastic entropy produc- 
tion Astot, the ensemble entropy production from Ref.— 
is recovered. 



D. Stochastic Gibbs relation 



(36) 



We can now relate the identified first law (f2"0")) to the 
usual equilibrium Gibbs relation. Inserting the definition 
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of the medium entropy (|30|) into ([20]) . we find 

AE = -TAs m - M« (42) 

In equilibrium and for quasistatic transitions, we have 
Astot = and therefore As = — As m which implies 



AE = TAs + y^fXjAr 



(43) 



Specialized to equilibrium, our approach thus yields 
along a single stochastic trajectory the usual thermody- 
namic relations valid for ensemble averages in equilib- 
rium. 



V. FLUCTUATION THEOREMS 

Fluctuation theorems provide exact results for driven 
systems beyond linear response. So far, they have been 
studied extensively for Langcvin dynamics and for master 
equations in general. Previous applications to chemical 
reaction networks^ 7 -^ have neither made the connection 
to an energy balance nor considered the entropy on the 
stochastic level with occasional exception o 31 ' 49 ' 50 ' 54 . The 
purpose of this section is to develop this concept in full 
generality. 



A. Integral fluctuation theorem 

The weight p[n(r),c(r)] for a particular trajectory 
n(r) with given initial state no under the concentration 
protocol c(t) is given by 



p[n(r),c(r)] = J[ exp - ^ 

i=l L r,p 

N, 

i=i 



dr<(n ; ,c(t)) 



with T[ (I — 1, . . . , TV;) being the jump times, where a 



reaction pi in direction r\ 



from state n, to state 



n, = n, ± v 



(45) 



occurs. For ease of notation, we have introduced the 
initial time tq ee 0, the final state n^- j+1 



final time tn 



R ee In 



t. We then define 

p (n ) •p[n(r),c(r)] 



iipf and the 



p(n t ,t) ■p[tl(t),c(t) 



(46) 



as the logarithm of the ratio of trajectory weights for 
a given trajectory n(r) with initial state no = n(0) 
and final state n t ee n(t) and the corresponding trajec- 
tory n(r) = n(t — r) under the time reversed protocol 
c(t) = c(i — t). Here, po( n o) is the initial distribution 



and p(n t ,t) is the solution of the master equation at time 
t. The integral fluctuation theorem 



<e-«) = 1 



(47) 



can then be proven by the following lines of 
identitie s 22 ' 27 ' 28 ' 31 using the normalization condition for 
the trajectory weight (j44|) and the final distribution 

1 = 5>(n t ,tMn(T),c(r)] 

n(r) 

- R ^p (n )p[n(r),c(r)} 



n(r) 



-fl[n(r)] 



n(r) 



Po(n )p[n(r),c(r)] 



-fl[n(r)] 



(48) 



Our thermodynamic approach allows to show that the 
abstract quantity R is exactly the total entropy produc- 
tion Astot- Inserting the trajectory weights lHH m t° ® 
yields 



R ee In 



hi 



p (n )p[n(T),c(T)] 
p(n t ,t)ij[n(r),c(r)] 



Po(n ) 



p(n t ,t) ' ^"^(n+d)' 
where c; ee c(t;). We then use (|3"Tj) to get 
p (n ) x > . /nr! 



(49) 



R = In 



p(n t ,i) 



5>(-^(n/ w ) 



As + As m = Astot- 



E 1 * 



<. ! (n; ,c/) 



wf! ri (nr,cj) 



(44) Here, we have used the shorthand notations 



and 



(£l/uj) n EE (Q/LJj)™ 3 



(50) 



(51) 



(52) 



The key point in this analysis is the cancellation of the 
degeneracy terms from the change in system entropy 
As with those from the entropy change of the medium 
As m in eq. (|50[) . Thus, the total entropy production 
Astot = R depends only on inital and final distribution 
Po( n o)> p( n tj*) and transition probabilities u/? _(n, c), 
see eq. (|49j) , which can, in principle, be extracted from 
experiments. Explicit knowledge of the normalizing vol- 
umes uj a and tUj defined only via the microscopic Hamil- 
tonian is not required. 
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The fluctuation theorem 



( e - As '°') - 



1 



(53) 



then constrains the probability distribution of the total 
entropy production As to t from which the second law 



(As to t) > 



(54) 



follows directly via the Jensen inequality (e x ) > . 
The relation (|53p holds for any initial state, any time- 
dependent protocol Cq,(t) and any length t of trajectories. 
Similar relations have been derived within a Hamiltonian 
dynamics in Ref^. 



B. Nonequilibrium steady states and detailed 
fluctuation theorem 

Likewise, for nonequilibrium steady states with con- 
stant c Q , the detailed fluctuation theorem 



p(-As t0 



-As tt 



XAstot) 



(55) 



can be shown exactly in the same fashion^. Here, 
p(Astot) is the probability density function for the to- 
tal entropy production. For the validity of this relation 
in the steady state for finite length t of trajectories, the 
inclusion of the stochastic entropy of the system As is 
crucial. In the long run, As m scales as t whereas As 
remains bounded. 



C. Generalized Jarzynski relation 



The Jarzynski relatio n 20 ' 21 



(56) 



expresses the free energy difference AF between two equi- 
librium states with a non-equilibrium average of the work 
TU[:e(t)] spent in this transition, which is a functional of 
the system trajectory x(t). In the form (|56"|) , it is valid 
only if the work is defined as 



i. Jo 



dE 



-A, 



(57) 



where the time-dependent protocols Afc(r) (equivalent 
to the c a (r) in our context) determine the transition 
rates via the time-dependent potential E(Xk{r), x). This 
definition of work is usually appropriate in a canoni- 
cal ensemble. We now derive a generalized version of 
the Jarzynski relation for our grandcanonical system. 
We consider transitions between equilibrium states and 
rewrite the quantity As to t using ([20]) and (|3U)) as 



As to 



As + PQ = As + f3W ( 



ch e 



(3 w chc 



J2w 



(3AE 



AJ (58) 



where we have used the definition of the grandcanonical 
potential 



A J = AE - TAs + Mj Ar 



(59) 



We now show that this change in the grandcanonical po- 
tential does not depend on initial (no) and final (n t ) state 
of the trajectory but only on the initial and final values 
of c. Using the grandcanonical equilibrium distribution 
and the definition of system entropy (|39p we get 



./ = E -Ts-^fijiij 

3 

= J2(Ej - Hj)n 3 + Tlnp cq (n, c) - ln.g n + lnW(c) 

3 

= J2( E 3 ~ MjK ~ J2( E J n i ~ + ln N(c) 

j j 
= ku\f(c) (60) 

and therefore A J = lnA/"(c t ) - lnA/"(c ). Using ([55]), the 
fluctuation theorem turns into the generalized Jarzynski 
relation 

^ e -/3[W chom -A(E jW n J )]^ = g -/3AJ (gl) 

Two remarks on the generality of this relation are appro- 
priate. 

(i) At constant c° q belonging to an equilibrium state, 
the exponent on the lhs of cq. (foTj) can be transformed 
to 



Wchom - A f ^ fljllj 

= J*dr (^2MC)n)j -A feton^j 



= 



(62) 



using the definitions of chemical work IU c hem ([2Tj) and 
chemical potential \ij (|2T)l . Thus, the result (|6"Tj) stays 
valid even if we do not await relaxation to the final equi- 
librium state. Such final relaxation implicitly assumed in 
the derivation of (|6"Tj) does not contribute to the exponent 
on the rhs. 

(ii) If we could clamp the initial and final number of 
molecules no and n t , we would get the original relation 



,-m\ - e -i 3AF , 



(63) 



D. Generalized Clausius inequality 



As mentioned above, there are two types of nonequilib- 
rium in chemical reaction networks : (i) a genuine NESS 
and (ii) any type of time-dependent c a (r) . It is useful to 
split 



As to 



In 



Po(np) 



£ln<^ 

l 



w p l ri (n+,c ; ) 



into a (c Q -dependent) part 
r ^ ln Po(n ) 



>(n t ,t) 



E 



Y + /3Q hk 
(64) 

(65) 



which accounts for non-adiabatic transitions and the so 
called housekeeping hea t 30 ' 49 ' 51 



iik 



E ln - 



(n ; ,c;)p s (n ; ,c ; ) 



(66) 



which is the permanently dissipated heat in a nonequilib- 
rium steady state for constant c a . Obviously, Qhk = if 
the detailed balance condition © holds for the station- 
ary state p s (n,c(r)) for any time r. For systems in a 
steady state, Y vanishes. For transitions between steady 
states, i.e., starting and ending in a stationary state, the 
quantity Y can be transformed to 



Y 



11 ^ 

Jo 



<9$ 

dc a 



(67) 



with $ = — hxp 8 . It has been show n 30 ' 49 , that the quan- 
tity Y obeys the integral fluctuation theorem (Hatano- 
Sasa-relation) 



(68) 



From (|68j) we get a generalized Clausius inequality even 
for transitions between nonequilibrium steady states 



(Y)>0 

or 

d (Qex) > - (As) 

with the excess heat 



(69) 



(70) 



0Q e 



As n 



hk 



7[ P s W,ci)g n + 



= Y - Ai 



(71) 

The relation ([70)) is a generalization of the classi- 
cal thermodynamic Clausius inequality Q > —AS to 
transitions between nonequilibrium steady states fit- 
ting in the phcnomenological framework of steady state 
thermodynamics 4 ^!. The fluctuation theorem 



1 



(72) 



for the housekeeping heat has recently been proven for 
stochastic processes ruled by a Fokkcr-Planck equation^. 
It can be shown that this relation holds also for chemical 
master equations^. 




FIG. 2: Model reaction network (|73|l with three freely evolv- 
ing species Xi, X2, X3 and two chemiostat species A and B. 



VI. AN EXAMPLE 

We illustrate our general results by considering the 
simple model reaction network (see Fig. [2J 



A- 



Xi X2 



k + 

X2 ■<— A3 



fe+ 

A3 ^- Xi 



B 



(73) 



where A and B are chemiostat species (corresponding to 
A\ and Ai in our previous notation). For simplicity, we 
assume that all normalizing volumes Wj (j = 1,2,3) in 
principle defined through a microscopic Hamiltonian, sec 
App. El arc equal. These normalizing volumes connect 



the bare rates fc+ _ to the transition rates 



.(n, c) 



through ([5]) and ([5]). The chemiostat species have chem- 
ical potentials (|15[) 



and 



HA = E A +Tln (c a uja) 



fJ-B = E B + T In (c b ojb) 



(74) 



(75) 



where Ea and Eb are single particle energies of A and 
B species, respectively, and u>a,ojb are appropriate nor- 
malizing volumes. The stochiometric matrix V = (v?) 



V = 




(P=l) 
(P = 2) 
(P = 3) 



(76) 



gives the number of transformed Xj molecules in one 
forward reaction of type p. Obviously, it is singular with 
rank V = 2, which corresponds to case II from Sect. 
Ill Bl We therefore have the constraint 



«2 



n 3 = N 



(77) 



on the accessible states of the phase space. Thus, the 
total number N of A-molcculcs is a conserved quantity. 
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A. Equilibrium distribution 

In order to determine the potential equilibrium distri- 
bution, we try to solve eqs. (CO 



-\h(caoja) + ln(fc_/fc + ) 

hi(k-/k + ) 
ln(c B w B ) + ln(fc_/fc+) 
(78) 




which have solutions only if 

ca/c b = [lob/ua) (fc_/fc+) 3 = 



cq 



B ■ 



(79) 



This is the expected condition for equilibrium since for 
any network cycle obeying detailed balance the product 
of forward and backward rates must be equal. We then 
get the solution 



ln(c A uJA) - ln(fc_/fc+) 


ln(/c_/fc + ) 



with arbitrary 7. The equilibrium distribution is the pro- 
jection of the Poisson distribution (fT0|) to the accessible 
state space, i. e. the states with total number N of 
X-moleculcs. This yields 







: 


H 








^(n,c) oc ft e -w 



+n 2 +n 3 



k+CAUJA 



hence the bare rates k + and k_ were equal. The ratio of 
bare rates k_/k + thus depends on the reservoir energies 
through (|83|) . We then get for the relationship between 
energies Ej and bare rates 







(!) 


'-( 







Eo — Ea 
E + ln(k-/k+) 
E + 2\n(k_/k+) 



(84) 



with arbitrary Eq. 

Having obtained single particle energies in an equilib- 
rium situation, the reaction system can now be driven 
out of equilibrium by changing the concentration of A 
or B molecules. Non-adiabatic transitions between equi- 
librium distributions are achieved by driving ca{t) = 
(lob/ua) (/s-//j+) 3 cb(t). Nonequilibrium steady states 
and transitions between them are obtained for any pro- 
tocol with c a (t) ^ (ujb/u>a) (A:_/fc + ) 3 c B (r). 



C. First law 

Whenever a single reaction takes place, the 
chemiostats associated with A and B molecules 
apply work to the system consisting of the Aj-molecules. 
Furthermore, the system of the A^-molecules dissipates 
heat into the surrounding aequous solution acting as 
a heat bath. The chemical work Wchem done by the 
chemiostats due to a single reaction is generally given 
by HU) as 



^chem 



r,p 



(85) 



B. Energy 

In principle, the energies Ej of one molecule Xj are 
defined as the thermodynamic average of the microscopic 
Hamiltonian, see App. [S] Without specifying such a 
Hamiltonian, for given bare rates the energies Ej 

of a single solute molecule Xj must obey (|2U1) 




/ln(fc_/*+) 
ln(fc-/*+) 
Vln(fc_/fc + ) 




which has a solution for the Ej only if 
E B -E A = 3ln 1 ^. 



(83) 



This is the condition for the reaction network to be ther- 
modynamically consistent as discussed in Sect. MI CI To 
simplify matters, we have set here and in the following 
the inverse temperature (3 = 1. 

The difference in energies Ea and Eb, combined with 
entropic forces, is the driving force of a potential nonequi- 
librium steady state. If the energies were equal, there was 
no reason for the reaction to run on average in one direc- 
tion for equilibrated concentrations caoja = cb^b and 



Due to energy conservation, the dissipated heat for one 
reaction step is given by (fT9)l as 



^EjAn^. 



(86) 



Having defined work, dissipated heat and internal energy, 
we can formulate the first law in the form (f^Uj) 



chem 



Q 



(87) 



Chemical work, dissipated heat and change of internal 
energy for one single forward reaction step of each type 
(p = 1, 2, 3) are shown in Tab. [J The expressions for the 
heat agree with the general result (f3"Tj) . The dissipated 
heat for a given sequence of reactions thus depends on the 
concentrations ca and cb- In equilibrium, the dissipated 
heat is zero for a complete cycle and therefore is zero on 
average. However, for a single reaction step the heat is 
not necessarily zero. 



D. Entropy production 

For an illustration of the entropy change, we calculate 
the total entropy production (|4"Tj) due to one occuring 
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Q 



AE 



HA 



-E 2 + Ei + [i A 
= in ( k + CA 



]n(k-/k+) + E A 



X 2 







-Ez + E2 



In 



X 3 Xi + B 



-E! 

= In 



+ Ez — /is 



-E A - 21n(fe_/fe+) 



TABLE I: Work, dissipated heat and change of internal energy 
for a single step in forward direction for the reaction network 
(ED). 



forward reaction 



X. 



fe_ 



B. 



(88) 



from state n = (ni, 712,^13) to state n + = (jii + 
1,712,713 — 1) at time r. The total entropy production 
Astot (EH) is the sum of the change in system entropy 



As = -lnp(n + ,r) + s°(n + ) + lnp(n , r) 
p(n~,r)n 3 
%(n+,r)(ni + 1) 

and the entropy change of the heat bath 

k 

As m — Q — In ' 



The total entropy production then is 



An") 
(89) 

(90) 



As tot = As + As m 
p(n~,r) 



= In 



In - 



In 



p(n+,r) (ni + l)fc_c B w B 
p(n~, -r)w+(n~, c(t)) 



p(n+, T)w^(n+, c(t)) 



i?. 



(91) 



This illustrates the cancellation of the degeneracy terms 
and the equivalence of the abstract quantity R ([^l 
and the thermodynamically motivated entropy produc- 
tion Astot SI])- The total entropy production thus could 
be extracted from experiments if the state probabilities 
and the transition probabilities of the stationary state 
could be measured. 



VII. CONCLUSION AND PERSPECTIVES 

In summary, for arbitrary chemical networks in a di- 
lute solution described by a chemical master equation, 
we have developed a consistent thermodynamic descrip- 
tion along a single stochastic trajectory by identifying 
work, internal energy and dissipated heat for each reac- 
tion step. A unique identification of internal energy and 
dissipated heat requires a microscopic Hamiltonian for 
the interaction between the reacting species and the sol- 
vent and cannot be extracted from the stochastic dynam- 
ics alone. However, we have shown that work and total 



entropy production do not depend on the microscopic de- 
tails and can be extracted from the dynamics once the 
transition rates in the chemical master equation are spec- 
ified or measured in an experiment. Entropy production 
involves both the change of the entropy of the medium 
given by the dissipated heat and genuine entropy of the 
network which requires a proper consideration of a de- 
generacy factor. The total entropy production fulfills an 
integral fluctuation theorem for arbitrary driving and a 
detailed fluctuation theorem in the steady state. 

This systematic approach relies on three crucial con- 
cepts. First, an energetic interpretation of stochas- 
tic master equation dynamics requires the comparison 
with an appropriate equilibrium state whose microscopic 
Hamiltonian for the interaction with the solvent must 
be known. By choosing the equilibrium concentrations 
of reservoir particles and thus adjusting the net flux to 
zero, a stationary distribution obeying detailed balance 
can be obtained for any thermodynamically consistent 
reaction network. Second, the persistence of mass ac- 
tion law kinetics even in non-equilibrium leads to the 
concept that energy differences are related to the ratio 
of bare rates in the same way as in equilibrium. Cru- 
cial for the persistence of this assumption is the time 
scale separation of diffusion and reaction time constants. 
For higher concentrations or long range interactions, a 
coupled system of reaction-diffusion equations would be 
necessary. The general concepts of stochastic thermo- 
dynamics would still apply; however, the expressions for 
energy, dissipated heat, work and entropy would get more 
complex. Third, we have decided to define the system as 
the collection of those species whose numbers we can fol- 
low in principle. The internal energy then is associated 
with these molecules only. This third assumption could 
be somewhat relaxed since one could always include one 
of the chemiostat species to the system species. The first 
two assumptions, however, are crucial if one wants to 
keep the energetic interpretation. Of course it is still 
possible to derive integral and detailed fluctuation theo- 
rems for a stochastic entropy production in any network 
obeying master equation dynamics 3 -^. The interpretation 
of the abstract quantity As m as an exchanged heat, how- 
ever, requires the first law energy balance and the correct 
identification of the degeneracy factor. 

Distributions of work or entropy production may be 
used to reconstruct thermodynamic equilibrium quanti- 
ties or to characterize the network, similar to the recon- 
struction of free energy landscapes for unzipping or un- 
folding transitions using the Jarzynski relatio n 32 i 35 ' 36 . A 
first step in this direction is the use of fluctuation the- 
orems to determine chemical driving forces in reaction 
cycles^. Of particular interest would also be the deriva- 
tion of the distributions of work and total entropy pro- 
duction for frequently used models of gene regulation net- 
works or signal processing networks, especially since it is 
known from a simple athermal model system that these 
distributions can show a quite rich structur o 59 i 60 . 

While there exists already a large amount of ex- 



12 



pcrimental work on fluctuation theorems for Langevin- 
systems, we are not aware of any attempts to probe fluc- 
tuation theorems for chemical reaction networks exper- 
imentally. Such experiments could contribute substan- 
tially to the ongoing effort to achieve a better under- 
standing of nonequilibrium systems. 
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APPENDIX A: NORMALIZING VOLUMES AND 
THERMODYNAMICS OF A DILUTE SOLUTION 

For a precise definition of the normalizing volumes 
uj a , ujj and the energies E a , Ej , we will now briefly re- 
view the standard thermodynamics of a dilute solution 
as found in textbooks, see e. g.— . For ease of notation, 
we treat the case of a dilute solution of N' equivalent 
solute molecules in a solvent with N molecules. This 
calculation can easily be generalized to the case of more 
than one species and applies both to the A a molecules 
and to the Xj molecules. We assume that the number 
of solvent molecules ./V is fixed and therefore use a semi- 
grandcanonical ensemble. In the dilute limit N' <C N, 
the Hamiltonian can then be split into the contributions 



H N > N '({r},{r'}) = 



with 



<({!■}) + <'({r'}) 

+H?- N '({r},{r'}) 



N' 



<'({r'}) ^E^«) 



(Al) 



(A2) 



and 



iV' 



H?> N ({r'},{r})=J2 H ^{ r })> ( A3 ) 

i=l 

where r[ denotes the position of the ith solute particle 
with i = 1,...,N' and {r} are the coordinates of the 
solvent molecules. One particle energies Ei (i = a,j) are 
defined through 



N,N' = 1 



(A4) 



Here, (.) denotes the canonical average with the appro- 
priate canonical Boltzmann- weight. 

Having defined one particle energies E a from a micro- 
scopic point of view through (|A4[) . one can transform the 
chemical potential @ into the form (fT5|) by choosing an 
appropriate normalizing volume 



ui a = e 



/co- 



CAS) 



where /3 = 1 /T is the inverse temperature. 

We now recall how to obtain the grandcanonical equi- 
librium distribution p e ^(N'). Integrating out the spatial 
degrees of freedom, we get 



yN,N' p f3fj,N' 
^can e 

V z N ' N 'pl 3 ^ N ' 



(A6) 



where Z^an is the canonical partition function which 
can be calculated as 



rN,N' 



N 



(A7) 



where fy(T,Q/N) and F(N) depend on the particular 
Hamiltonian and A' is the thermal de Broglie wavelength 
of a solute particle^. The solute molecules thus be- 
have like an ideal gas in the effective (constant) potential 
V(T,n/N)/[3. Inserting flST]) into {SjBJ yields 

N' 



1 



A" 



-fi(E-n)N' 



(A8) 



Thus, Eq. jXl]) has the form of cq. ((28]) with 

cu = \' 3 e- pE /y. (A9) 

A simple calculation shows that eqs. (jTSJ) and (|A5[) are 
consistent with eq. (|A9[) . 



APPENDIX B: ANALYS IS OF CASE II FROM 
SECT. [TIB] 



If there exists more than one solution of (fT2"|) , the rank 
of the stochiometric matrix V is smaller than the number 
Nj of species Xj . This is usually the case if there are less 
reactions p than species Xj, i.e. N p < Nj. We then 
rankY (generically : N p = Nj — N p ) 



have N p = Nj 
constraints 



E 



(Bl) 



fi = (l,...,N p ) for the accessible states with L = 
(L , . . . , depending on the initial distribution. The 
matrix K has rows which span the nullspace of the sto- 
chiometric matrix V 



(B2) 



The general solution of (jT2j) then is the sum of a special 
solution e°j and a nullspace vector Kx 



E A 7 a 



(B3) 
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for an arbitrary vector x 6 R M . The state space {n} 
is then separated into an infinite number of subspaces 
Xl = { n \K jUj = L M } C {n}. If all possible initial states 
lie in the same subspace xl, the equilibrium distribution 
is the projection of the Poissonian (fT4]) to the subspace 
XL with the Cj given by any solution of (fl2|) . It does not 
matter which solution of (JT3J) is chosen since 



Projections comprise normalization and therefore we get 
the same distribution for different solutions of (fT2)) . 



P eq (n,c) = 



rijl 



i J 
^"Q (^M)" J c -( e , n ,)-L, 



(B4) 
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